Cost function for low-dimensional manifold topology assessment

In reduced-order modeling, complex systems that exhibit high state-space dimensionality are described and evolved using a small number of parameters. These parameters can be obtained in a data-driven way, where a high-dimensional dataset is projected onto a lower-dimensional basis. A complex system is then restricted to states on a low-dimensional manifold where it can be efficiently modeled. While this approach brings computational benefits, obtaining a good quality of the manifold topology becomes a crucial aspect when models, such as nonlinear regression, are built on top of the manifold. Here, we present a quantitative metric for characterizing manifold topologies. Our metric pays attention to non-uniqueness and spatial gradients in physical quantities of interest, and can be applied to manifolds of arbitrary dimensionality. Using the metric as a cost function in optimization algorithms, we show that optimized low-dimensional projections can be found. We delineate a few applications of the cost function to datasets representing argon plasma, reacting flows and atmospheric pollutant dispersion. We demonstrate how the cost function can assess various dimensionality reduction and manifold learning techniques as well as data preprocessing strategies in their capacity to yield quality low-dimensional projections. We show that improved manifold topologies can facilitate building nonlinear regression models.

In the era of big data, numerous science and engineering disciplines use dimensionality reduction to obtain lower-dimensional representations of complex physical systems with many degrees of freedom [1][2][3][4][5][6][7][8] . Large data coming from system measurements or simulations is frequently the starting point of reduced-order modeling. These high-dimensional datasets, ubiquitous in areas such as plasma physics, chemically reacting flows, neuroscience, genomics and transcriptomics, electrochemistry or atmospheric physics, often exhibit strongly attracting low-dimensional manifolds [9][10][11][12][13][14][15][16][17][18][19][20] . Describing the system evolution on those manifolds alone can thus be a viable modeling strategy [21][22][23][24] . After projecting the original variables onto a lower-dimensional basis, system dynamics can be tracked on a lower-dimensional manifold, embedded in the original state-space. This approach provides substantial reduction to the number of parameters needed to visualize, describe and predict complex systems. To date, linear and nonlinear dimensionality reduction techniques have been used to find lower-dimensional spaces to represent multivariate datasets and build reduced-order models in those spaces.
Some topological properties of low-dimensional data representations can make reduced-order modeling difficult. A particularly undesired behavior are overlapping states on a manifold which can result in non-uniqueness in dependent variable values. For instance, with the manifold parameters used as regressors in nonlinear regression tasks, any ambiguity in dependent variable values can hinder successful modeling. Another characteristic of a problematic manifold are large gradients in the dependent variable values. These can appear when observations on manifolds in important regions are compressed with respect to other, less important regions. If dependent variable values change rapidly over such compressed regions, features of small sizes are formed that can pose modeling difficulty 25 . Maintaining moderate gradients on manifolds is thus a desired characteristic.
In this paper, we are motivated by the emerging discussion on the need to characterize the quality of lowdimensional manifolds [26][27][28][29][30] . Quantitative tools are needed in areas where researchers tune the hyper-parameters of dimensionality reduction or manifold learning techniques to obtain improved manifolds of particular types of data [31][32][33][34][35] . When manifold learning is used for efficient data visualization, good quality manifolds can help uncover important multivariate relationships in complex datasets such as genomics or transcriptomics data 34,[36][37][38] . In artificial intelligence, there is a need to determine the quality of manifold representations in neural networks, where the recent work on predictive learning demonstrated how those representations can vary throughout the learning stages 19 . Detecting intersection between several manifolds can be of interest in learning object manifolds

Results
Cost function formulation. We base our discussion on the observation that regions of poor manifold topology (such as non-uniqueness) will only affect a dependent variable if there is variation in the variable's values over those regions (see Fig. 1). With this premise, our cost function is derived from variance in dependent variable values happening across different length scales on a manifold. As recently proposed in Ref. 30 , we compute the normalized variance of a dependent variable, N (σ ) , for selected length scales, σ , on a manifold. We then analyze the derivative of N (σ ) with respect to σ , denoted as D (σ ) . This derivative captures the information about how the normalized variance changes as the manifold is scanned at varying length scales given by σ ∈ �σ min , σ max � . The detailed mathematical description of N (σ ) and D (σ ) can be found in the "Methods" section and in Ref. 30 . Given q manifold parameters, η η η = [η 1 , η 2 , . . . , η q ] , the resulting normalized variance derivative, D (σ ) , can be computed for each dependent variable in a relevant set, φ φ φ = [φ 1 , φ 2 , . . . , φ m ] . Figure 2a shows a visual example of how the normalized variance derivative assesses information content at various manifold length scales for one dependent variable, φ . The length scales at which peaks occur in a variable's D (σ ) profile Various 2D projections of a 3D synthetic dataset can be formed by looking at the dataset at various angles to the z-axis and collapsing all observations onto a plane of sight. We demonstrate two example projections that can be formed: a top-down projection and a projection resulting from looking at an angle to the z-axis. In the new projected coordinates, [η 1 , η 2 ] , the top-down projection is unique everywhere, while the projection at an angle introduces regions of non-uniqueness with overlapping observations. www.nature.com/scientificreports/ indicate feature sizes. Each scale at which a dependent variable shows variation over the manifold, including variation due to non-uniqueness, will create its own imprint in D (σ ) . The largest spatial scale at which the dependent variable exhibits variation (the largest feature size) on a manifold is reflected in the rightmost peak location in the D (σ ) profile, and we define it as σ peak (see Fig. 2a). In the "Methods" section, we present a more detailed discussion of how σ peak is obtained. Any smaller feature sizes appear as additional peaks for σ < σ peak . The 2D projection seen in Fig. 2a exhibits severe non-uniqueness which manifests itself in variance occurring at very small scales, σ ≪ σ peak . Two particularly appealing characteristics of the normalized variance derivative are taken into account for the design of our cost function. First, the locations of peaks in the D (σ ) curve convey feature sizes on manifolds. With an appropriate penalty with respect to the largest feature size, σ peak , we favor manifolds that maintain large feature sizes, as those should facilitate modeling. Second, with multiple scales of variation present on manifolds, the area under the D (σ ) curve will increase due to the additional peak(s) for σ < σ peak . The cost function proposed herein, L = L(η η η, φ φ φ) , is computed from a penalized area under the D (σ ) curve(s). By integrating the penalized D (σ ) curve(s), we sum over the effects that multiple scales of variation have on the D (σ ) profile. A visualization of the penalty function and its effect on the D (σ ) curve for the manifold introduced in Fig. 2a is shown in Fig. 2b. We see how the area under D (σ ) weighted by the penalty function is amplified at scales far from σ peak , therefore penalizing the variance occurring at smaller scales more heavily. Furthermore, the cost function gives a single number representing the parameterization "cost" for a given manifold. A larger cost indicates a worse manifold topology and a lower cost indicates an improved manifold topology in terms of modeling. If the cost is computed for m dependent variables, a norm over all L i costs can be taken to yield a single cost value for the entire manifold: L = ||L i || , ∀i ∈ [1, 2, . . . , m] . Figure 2c demonstrates costs computed for a single dependent variable, φ , over two different 2D projections. The first projection is the same as analyzed in Fig. 2a. The second projection has an improved topology and non-uniqueness is significantly reduced. The corresponding D (σ ) curve for the improved projection exhibits a single dominant peak which indicates a single scale of variation in φ values over the entire manifold. We report the cost, L = L φ , for each projection, with the cost for the projection with non-uniqueness being greater. The mathematical description of the proposed cost function and additional details are provided in the "Methods" section.
There are three key advantages of our proposed cost function. First, manifolds obtained using any technique can be assessed. This includes any ad hoc selected manifold parameters or empirical manifolds obtained directly from training data using dimensionality reduction or manifold learning techniques. Second, manifolds of any dimensionality can be assessed. Finally, manifolds can be assessed with respect to an arbitrary set of m relevant dependent variables, φ φ φ . Dependent variables that we are most interested in accurately modeling can be selected, which can include any of the original state variables and any functions of them.
Cost function response to feature size and non-uniqueness. To demonstrate the behavior of the proposed cost function to an increasing feature size, we generate bivariate functions on a uniform square xy grid, centered around (0, 0). The dependent variable, φ , is calculated from a multivariate Gaussian normal distribution, φ = exp(−(x 2 + y 2 )/(2s 2 )) , where the standard deviation, s, is gradually increased to imitate an increasing feature size. The smallest feature is selected with s = 0.05 and the largest feature with s = 0.6 . Figure 3a shows these functions corresponding to ten Gaussians with increasing s. Below each Gaussian, we plot . The location of the rightmost peak, σ peak , denotes the largest feature size on a manifold. Variance occurring at scales σ ≪ σ peak is a strong indicator of non-uniqueness. (b) The cost function, L , is computed by integrating the penalized D (σ ) curve. The introduced penalty amplifies the area under D (σ ) occurring at σ < σ peak and especially at σ ≪ σ peak . The contributions from all scales of variation are thus summed up and embedded in the cost value. (c) Assessing example 2D projections with the proposed cost function. The first projection exhibits non-uniqueness, with regions where smallest and highest values of the dependent variable, φ , overlap each other. The non-uniqueness is significantly reduced on the improved projection and the corresponding D (σ ) curve exhibits a single rise indicating a single feature size over the whole manifold. We show the associated costs for the two projections, with the cost for the projection with non-uniqueness being greater.
the corresponding cost, L = L φ , and we observe a decreasing trend in L with increasing feature size. This is a desired behavior, since larger features should facilitate modeling. A very small feature, like the one obtained from s = 0.05 , introduces relatively steep gradients, which can be more challenging to model. The response of the D (σ ) curves to an increasing feature size is shown in Fig. 3b. The location of σ peak gradually shifts to the right with increasing feature size and the area under the D (σ ) curves decreases. The decrease in cost with an increasing feature size seen in Fig. 3a is thus a result of rewarding the increasing σ peak location and the slightly decreasing area under the D (σ ) curve. The latter is the effect of a decreasing gradient of φ with increasing feature size, leaving less variance in φ present at scales away from σ peak . To further test the response of the cost function to multiple feature sizes, we generate another set of functions, such that the dependent variable, φ , is now computed from a superposition of sine functions with various frequencies. Multiple feature sizes are generated from the general formula φ = n k=1 sin(2 k x) for n = 1, 2, . . . , 5 . Thus, for n = 1 , the function has only one feature size. In Fig. 3c, we observe increasing cost values with each added feature size. The reason for this increase becomes clear once we look at the corresponding D (σ ) curves in Fig. 3d, where each new feature generates an additional rise in D (σ ) at length scales smaller than σ peak .
Next, we perform a similar test using functions which introduce increasing levels of non-uniqueness. The dependent variable, φ , is now calculated as a linear function of the independent variable, φ = x , being the x-axis. This results in a constant gradient along the x-axis. We introduce non-uniqueness in φ by adding overlapping observations whose φ values are set to zero. The number of observations added to increase the non-uniqueness is measured with the overlap depth, d, and is increased from no overlapping observations ( d = 0 observations) to a maximum number of overlapping observations ( d = 90 observations). Figure 3e presents ten such functions with an increasing depth of non-uniqueness. Here, we observe an increasing trend in L as we increase the non-uniqueness depth. The multiple scales of variation introduced from non-uniqueness lead to an increased area under the D (σ ) curves, which is penalized by the cost function. In order to observe how severe the increase in area is under the D (σ ) curves associated with manifold non-uniqueness, Fig. 3f shows the D (σ ) curves www.nature.com/scientificreports/ corresponding to the ten functions with non-uniqueness. The first interesting observation is that the location of σ peak for all D (σ ) curves is almost identical. This is understandable, since the size of the main feature (the constant gradient) stays the same. The only reason for the increased area under the D (σ ) curves is the appearance of additional peaks at length scales σ ≪ σ peak . Only in the case d = 0 (no overlap) the D (σ ) curve exhibits a single rise. As long as any overlapping observations are introduced, an additional area shows up under the D (σ ) curve at σ ≪ σ peak . This additional area is most pronounced for the largest non-uniqueness depth, d = 90. Finally, it is instructive to discuss the scales, σ , at which we detect overlap in Fig. 3f. The highest rise in D (σ ) linked to non-uniqueness is happening at σ ≈ 5 × 10 −4 for all cases d > 0 . This exact value can be linked to the sample density for the ten toy functions with increasing overlap depth. For a normalized manifold, where x ∈ [0, 1] , the distance along the x-axis between data points in the unique region is σ ≈ 10 −3 . The overlapping observations are located in-between every other unique observation, such that the distance along the x-axis between any observation from the unique region and its nearest overlapping observation is of the order of 10 −4 . Thus, as the manifold is scanned with varying length scales, σ , once σ ≈ 10 −4 , the variation in a dependent variable values is captured between a single point from the unique region ( φ → 1 ) and a single point from the overlapping region ( φ = 0 ). This creates a sudden increase in the captured variance and shows up as the peak at σ ≈ 5 · 10 −4 in the D (σ ) curves in Fig. 3f. A more in-depth discussion of how scales of variation can be linked to data density on a manifold grid can be found in Ref. 30 .
Assessing data preprocessing strategies. We now turn our attention to multivariate datasets and their lower-dimensional projections. Each dataset considered next is represented by a matrix X ∈ R N×Q , where N is the number of observations and Q is the number of state variables. Since most often N ≫ Q , Q determines data dimensionality. We consider three disciplines that deal with datasets that are notoriously high dimensional: plasma physics, reacting flows and atmospheric physics. For example, the state-space of a reacting flow is described by the temperature, pressure and chemical composition. There, the high state-space dimensionality originates from many chemical species involved that can easily reach the order of hundreds. In this and the following sections, we demonstrate a few practical applications of the proposed cost function using those datasets.
Reduced-order models leverage the fact that multivariate datasets can often be successfully re-parameterized by a reduced set of low-dimensional parameters. To date, numerous techniques have been employed for dimensionality reduction of multivariate data. Those include linear techniques such as principal component analysis (PCA), independent component analysis (ICA), distance metric learning (DML) or linear discriminant analysis (LDA) and nonlinear techniques such as kernel PCA (KPCA), Isomap 57, locally linear embedding (LLE) and its variants 58 or autoencoders 59 . For the purpose of data visualization, t-distributed stochastic neighbor embedding (t-SNE) 55 and uniform manifold approximation and projection (UMAP) 60 are gaining popularity in various research disciplines 7,36-38,61-64 . A short summary of dimensionality reduction and manifold learning techniques explored in this work can be found in the "Methods" section.
Prior to applying dimensionality reduction, training datasets are often preprocessed. The most straightforward strategy is data normalization by centering and scaling each state variable. Other preprocessing approaches involve data sampling to mitigate imbalance in observation density, or feature selection. The effect of data preprocessing alone can have a large impact on the resulting low-dimensional manifold topology constructed from such data 53,65 . To demonstrate how significant those changes can be, Fig. 4a shows 2D PCA projections of a dataset describing argon plasma (N = 100,700, Q=36) 66,67 , where the state-space is spanned by the following variables: temperature of heavy species, T h , temperature of electrons, T e , and 34 species mass fractions that contain 31 electronic states of argon, two levels of ionized argon, Ar and Ar + , and electrons. The projections are generated with various scaling techniques applied on the dataset (see the "Methods" section for the summary of all scaling techniques used in this work). In the top row of Fig. 4a, the projections are colored by T e and in the bottom row by the electron mass fraction, Y e . Our visual analysis shows that preprocessing can affect manifold topologies significantly. Various topologies are expected to perform differently in reduced-order modeling due to changes in feature sizes or the level of non-uniqueness.
We quantitatively assess those changes and rank the preprocessing strategies in their capacity to generate quality manifold topologies. In Fig. 4b, we show costs, L , for the various 2D PCA projections visualized above, and for comparison, for the analogous 3D PCA projections (although not visualized). Here, L is computed as the L 1 -norm over the individual costs for the relevant dependent variables: T h , T e , Y e and Y Ar . Thus, in Fig. 4a, the projections are colored by the two target dependent variables to allow us for making a visual connection between how the projection topologies look like versus how the cost function assesses a given projection. We note that VAST scaling results in the lowest cost for 2D projections. Looking at the corresponding visualization in Fig. 4a, we can see that T e and Y e values change relatively smoothly across the manifold, even though the VAST projection introduces the "spike" region, where many data observations are compressed. For comparison, projections resulting from the remaining scalings introduce relatively steep gradients and overlap for a few of them, in some region of the projection, either for the T e or the Y e variable. Finally, we note that costs drop across all scaling strategies explored when the PCA projection dimensionality is increased to 3D.
We additionally explore feature selection (otherwise known as variable subset selection) as another strategy in the data preprocessing pipeline. Feature selection involves finding a meaningful subset, X S ∈ R N×S , of the original Q state variables ( S < Q ), and using this subset in data science or machine learning algorithms, instead of the full set of state variables 68 . This time, we use a reacting flow dataset describing combustion of syngas in air. In Fig. 4c, we show 3D PCA projections of an 11-dimensional dataset (N = 14,550, Q = 11), where the eleven dimensions are spanned by temperature, T, and ten mass fractions of chemical species, denoted Y i for species i. The different projections seen in Fig. 4c  www.nature.com/scientificreports/ dataset (see "Methods"). The bottom row shows 3D projections resulting from different scalings combined with feature selection applied to the original training data. All PCA projections are colored by the temperature. For the purpose of this demonstration, we use a newly developed feature selection algorithm that uses the cost function to guide the optimal selection of the variables from the original training data 65 . By minimizing the cost of the resulting data projection, we optimize the feature selection process from the point of view of manifold topology. The detailed description of the feature selection algorithm can be found in the "Methods" section. In Fig. 4d we show costs, L , tied to the two preprocessing strategies explored: just scaling the data (circles) versus scaling with feature selection (triangles). Here, L is computed as the L 1 -norm over the individual costs for the relevant dependent variables: The dependent variables were selected manually as the most important candidates in modeling. Black markers show costs corresponding to the 3D PCA projections visualized in Fig. 4c. For comparison, with gray markers we show costs for analogous 2D PCA projections. Three of these 2D projections corresponding to scaling only are visualized in Fig. 4e-the two worst projections (largest L ) and the best projection (smallest L ). We note that the worst 2D projection corresponding to �−1, 1� scaling is severely folded over itself which is likely the reason for the high cost. The Level projection also introduces significant non-uniqueness, with the low-temperature regions being represented over narrow geometry on a manifold. The best 2D projection corresponds to VAST scaling. As seen in Fig. 4e, this projection is much more unique compared to �−1, 1� or Level scaling 2D projections. However, there is a visible "twist" on the VAST 2D projection, which is untangled when a 3D projection is considered instead. This is a possible factor for a decreased cost between the 2D and the 3D PCA projection. The lowest costs for 3D projections happened for VAST (for scaling only) and Auto (for scaling with feature selection). The corresponding visualized projections are marked with thicker axes in Fig. 4c. These two projections are characterized by relatively well-spaced feature sizes and reduced non-uniqueness. Finally, similarly as we have observed for the argon plasma data, costs drop across all scaling strategies when the PCA projection dimensionality is increased from 2D to 3D. We also compare costs corresponding to the 3D projections visualized above (black markers) with the analogous costs of 2D projections (gray markers). The optimal manifold topology corresponding to the lowest L for each preprocessing strategy is highlighted with thicker axes in (c). (e) Visualization of 2D projections corresponding to three selected cases of only scaling the original data: �−1, 1� and Level scaling (corresponding to the two highest L ) and VAST scaling (corresponding to the lowest L). www.nature.com/scientificreports/ Our analysis from Fig. 4 reveals that appropriate data preprocessing can help in generating a better quality low-dimensional manifold. In particular, optimized feature selection can be beneficial as it decreases costs over the relevant dependent variables with respect to costs associated with only scaling the data. With the many possible data preprocessing techniques encountered in the data science community, the cost function allows for quantitative rankings without the need to analyze manifold topologies manually. This, in turn, can help fine-tune data preprocessing to a given dataset and a desired manifold dimensionality.

Scientific Reports
Detecting large gradients on manifolds. Using the proposed cost function, manifolds can be assessed with respect to different relevant dependent variables separately. If there is a good coverage of regions where the selected variables vary over the entire manifold, problematic regions on manifolds can be detected with greater reliability. Figure 5a demonstrates a 3D manifold embedded in a 9-dimensional state-space of a reacting flow dataset describing combustion of hydrogen in air (N = 13,468, Q = 9). The topology of this manifold is curved such that two regions where fuel and oxidizer originate from are brought closely together. This region is better visualized in the zoomed-in dashed box. When the manifold is colored by the hydrogen (fuel) mass fraction in Fig. 5a, we see a step-change in the φ = Y H 2 values over the considered region. This potentially undesired behavior can be detected with the cost function as it increases the area under the D (σ ) curve at length scales smaller than σ peak . The cost associated with the hydrogen mass fraction is L H 2 = 1.8 . For comparison, the temperature variable does not experience large variation within the considered region (Fig. 5b) and the cost associated with temperature is lower ( L T = 1.0 ). We observe a single peak in the D (σ ) profile corresponding to the temperature variable.
Regions where opposing physical phenomena meet over small length scales on a manifold (in this case fuel and oxidizer streams) can prove to be difficult to model accurately. For instance, some regression models are known to struggle in the presence of large gradients in dependent variable values 25 . Such problematic regions of compressed observations on manifolds can be detected by analyzing costs across several dependent variables. If there exists an important dependent variable that varies across the problematic region (such as the Y H 2 variable in Fig. 5a), it can help discard the problematic manifold topology. We note that the choice of the relevant dependent variables is important in assessing whether a given manifold topology is appropriate or not from the modeling perspective. Some variables, such as the temperature variable alone in the example shown in Fig. 5b, might not be effective at exposing regions of non-uniqueness in data projections, that can affect other relevant variables. We also note that the dependent variables for which the cost is computed can be arbitrary and do not have to be selected from the set of the original state variables as we have done here. They can equally be functions of the state variables. In reacting flow applications, these can for instance be the production rates of chemical species, related to the state-space through nonlinear Arrhenius expressions.
Manifold assessment across dimensionality. We now demonstrate how the cost function applies to manifolds of arbitrary dimensionality. This can be useful in determining an appropriate dimensionality for rep- www.nature.com/scientificreports/ resenting the variables of interest when using a reduction technique. The demonstration is first done using an experimental combustion dataset known as the Sandia flame D dataset 69 . This dataset contains approximately N = 57,000 observations each for temperature and various species mass fractions (Q = 10) over six different heights in a methane and air piloted jet flame. The experimental dataset also presents the opportunity to demonstrate the cost function behavior on manifolds containing noise. Figure 6a shows an example 2D PCA projection with Max scaling of the Sandia flame D data colored by temperature. This manifold is less structured than those seen in Figs. 4 and 5 for the numerically simulated data. With noisy data, such as seen from experiment, it can be difficult to visually assess the quality of manifold topologies, especially in higher dimensions. The proposed cost function should adequately assess such topologies since the D (σ ) curves remain smooth even for noisy data 30 . Figure 6b demonstrates how the cost function behaves with increasing dimensionality of a projection using PCA, up to the original dimensionality of the system. The two opaque curves correspond to the scalings that resulted in the highest and lowest costs for the projected experimental data. Also included is the cost for the original ten-dimensional manifold. For the most part, we see a decrease in the cost with increasing dimensionality. Increasing the dimensionality further after four dimensions appears to have negligible effects on the feature sizes of the optimizing variables since the cost function flattens out. At this point, the computational cost of adding dimensions to a model most likely outweighs any benefit seen from slight increases in feature sizes. Any differences between the various scalings also become less noticeable as dimensionality increases. We also see that rotating the manifold along new principal component axes, even without reducing dimensionality, can slightly reduce the cost compared to the original manifold. This is due to variance being measured along new coordinates. For comparison, in Fig. 6b we also include costs with increasing dimensionality for the Sandia flame F dataset (transparent curves), which is another experimental case that exhibits stronger effects of flow turbulence than flame D. As a result, flame F experiences conditions with local flame extinction, absent in the flame D dataset. Thus, flame F includes many more states of methane combustion than flame D and the manifold corresponding to flame F covers a wider range of attainable thermo-chemical states. For instance, for a specific stoichiometric condition, the flame can be either in an ignited or in an extinguished state. We thus expect an increased required manifold dimensionality for flame F to achieve the same parameterization quality as compared to flame D. This is indeed the case seen in Fig. 6b, where costs are generally higher for flame F than for flame D. The proposed cost function could be used in this way to determine an appropriate minimum dimensionality for the projection of data, with the goal of facilitating modeling of the optimizing variables 70,71 . Manifold assessment across various dimensionality reduction and manifold learning techniques. We now revisit a recent PCA-based reduction of a high-fidelity dataset obtained from delayed detached eddy simulation in atmospheric physics 46 . This data simulates the Cedval A1-5 wind tunnel measurements of atmospheric dispersion of a pollutant in the vicinity of a rectangular building 72 . The numerical data that we use here represents a planar slice through the computational domain, perpendicular to the downstream wind direction. The dataset has N = 18,540, Q = 16 and the original state-space is composed of variables common to atmospheric physics applications, such as the three velocity components, pressure, Reynolds number, or the turbulent viscosity, kinetic energy and dissipation rate. All parameters in the data have been averaged over a 3.5 s period. We use an outlier detection algorithm to remove outliers from the data (see "Methods"). The aim of the original study 46 was to predict the turbulent Schmidt number, Sc t , from the PCA-derived manifold parameters. This demonstrates an application where regression is the end goal for obtaining a data-driven correlation of a physical quantity with the manifold parameters, η η η . We thus seek such regression function f that Sc t ≈ f (η η η) . While in the original work PCA was used to reduce the state-space dimensionality, here we benchmark three additional techniques: UMAP, Isomap and t-SNE. Following a similar methodology as demonstrated in Fig. 4, we first searched for the best scaling option for each reduction technique. In addition, we introduce three novel scaling methods inspired by the variable stability (VAST) scaling 49 . We denote them S1, S2 and S3 (see Table 1). www.nature.com/scientificreports/ S1 scaling is an extension of VAST, where the effect of non-normality is considered by multiplying the standard deviation by the data kurtosis. S2 and S3 are variations of S1 obtained by replacing the mean value in the coefficient of variation by the maximum and the range of each variable respectively. In Fig. 7a, we compare qualities of 2D projections, corresponding to the best scaling, both visually and using our cost function. We report costs for the φ = Sc t only, since that was the modeled dependent variable. An appropriate scaling allowed us, in  for ANN prediction of Sc t across all four projection techniques (blue diamonds) using 3D data projections. The minimum L and the minimum MAE for each technique is marked with a shaded outline. Scalings that generally exhibit lowest costs (VAST, S1, S2, S3) also result in the smallest MAE. (d) Example 3D PCA projections resulting from applying two scaling options to the original data: S1 and Pareto scaling. For this dataset, S1 scaling allowed for generating a manifold where the dependent variable of interest, Sc t , has a smooth gradient across one of the manifold dimensions. Pareto scaling, exhibiting high cost, collapses most of the data observations onto a planar structure. (e) The D (σ ) curves corresponding to the 3D PCA projections with S1 and Pareto scaling. (f-i) Scatter plots of L versus MAE from ANN and kernel regression predictions of Sc t . We show predictions based on 600 different 2D PCA projections (f, g), and based on 600 different 2D t-SNE projections (h, i) as the independent manifold parameters. We test a few selected scaling techniques applied to the atmospheric dispersion data; legend applies to all figures (f-i). www.nature.com/scientificreports/ many cases, to find significantly better projections. For instance, the worst 2D PCA projection is associated with L = 5.2 , while the best (visualized) with L = 1.2 . Figure 7b shows the D (σ ) curves corresponding to the 2D projections visualized in Fig. 7a and to φ = Sc t . The D (σ ) curves are mostly composed of a single rise, suggesting that the non-uniqueness in the 2D projections has been remedied by selecting an appropriate data scaling.
Improved manifold topologies yield more accurate regression. Nonlinear regression can be used in combination with dimensionality reduction to provide mapping functions, f, between the manifold parameters and physical quantities of interest. Regression can help discover robust relationships between manifold parameters and physical quantities that can then be injected in computational models and simulations 18,44,46,73 . Techniques, such as artificial neural networks (ANNs), Gaussian process regression (GPR) or kernel regression are commonly used in this context. Below, we continue the example of the atmospheric dispersion dataset and demonstrate improvements in nonlinear regression performance when parameters of an improved manifold topology are used as regressors. We first train an ANN model to predict Sc t based on 3D PCA, UMAP, Isomap and t-SNE projections resulting from different data scaling options. The details on the ANN model used here are provided in the "Methods" section. The scalings ranking for these 3D projections, along with the mean absolute errors (MAE) for Sc t predictions using ANN are reported in Fig. 7c. We note that generally, the scalings which resulted in lowest costs (VAST, S1, S2 and S3) across all four reduction techniques exhibit lowest MAE. The minimum L and MAE for each technique is marked with a shaded outline. With the exception of t-SNE, the minimum L and the minimum MAE happened for the same scaling option.
The reason for a good ANN performance on those four scalings can be further understood by visualizing example projections. In Fig. 7d, we show 3D PCA projections of the atmospheric dispersion data, corresponding to the best scaling option (S1 with L = 1.3 ) and the worst scaling option (Pareto with L = 4.2 ) as per the ranking shown in Fig. 7c. The quality of projections changes visibly with the change of data scaling. For the S1 scaling projection, we observe a clear gradient of Sc t throughout the manifold, while for the Pareto scaling case, the 3D projection introduced significant scatter in Sc t values, with many observations compressed onto a nearly planar structure in the reduced 3D space. The latter projection thus introduces significant non-uniqueness in Sc t values and is more difficult to accurately regress over. The difference in costs for the two 3D PCA projections from Fig. 7d can be further understood by looking at the comparison of D (σ ) curves in Fig. 7e. Pareto scaling case creates much more variance present at smaller length scales which is consistent with our visual inspection of the projection.
More generally, we observe correlation between the verdict given by L and the nonlinear regression performance. In Fig. 7f-i, we show scatter plots of L versus MAE for varying manifold topologies obtained from the atmospheric dispersion dataset. Here, we focus on PCA and t-SNE as the two projection techniques. Each scatter plot takes into account six selected scaling techniques: Auto, Pareto, VAST, 0, 1 , Level and S1. For each scaling technique, we generate 100 distinct 2D projections through creating random variable subsets (feature selections) of the full dataset, before applying PCA or t-SNE. This allows us to have sufficiently many distinct manifold topologies (600 for each scatter plot) for a trend to emerge in Fig. 7f-i. In addition to ANN regression, we apply kernel regression of Sc t to observe if the correlation between L and MAE is still present for a different nonlinear regression technique. The details on the kernel regression model used here are provided in the "Methods" section. Although the trends of L vs. MAE are different between ANN and kernel regression, the correlation measured with Spearman coefficient is high in each case. For ANN regression we observe 96% Spearman correlation for 2D PCA and t-SNE projections. The correlation is lower for kernel regression as compared to ANN regression, with 90% and 93% Spearman coefficient for PCA and t-SNE respectively. In Fig. 7f-i, manifolds corresponding to Level scaling are clustered in the region of high L and high MAE. This result is consistent with Fig. 7c, where Level scaling is seen to yield high costs and high regression errors. Conversely, the region of smallest L and smallest MAE in Fig. 7f-i is occupied by topologies resulting from S1 scaling, although some random subsets can result in poorer manifold topologies even with S1 scaling. We observe similar correlation trends for 3D PCA and t-SNE projections; these are shown in the Supplementary material. With the regression hyper-parameters kept constant throughout this exercise, the results presented in Fig. 7f-i suggest that better regression performance can be achieved when manifold topologies are improved. The cost function can thus be used to find optimized regressors for nonlinear regression techniques such as ANN or kernel regression. Other regression techniques are expected to be similarly affected.
Detecting overlap between classes in categorical data. So far, we have shown examples where the cost was computed using continuous dependent variables. Here, we briefly explore the cost function application to categorical data. The dependent variable, φ , can now be formed from numerical values of the class labels. Figure 8a shows in an illustrative way the capability of the D (σ ) metric to detect overlap between classes on a projection. The overlap between the two clouds of points, each representing one class, is reflected in an additional peak seen in the D (σ ) profile. When the clouds become sufficiently separated in space, D (σ ) exhibits a single rise. This behavior is then translated in the cost value, with the cost for the case with overlapping classes being greater. Next, we generate 2D projections of the MNIST handwritten digits dataset 74 using PCA and t-SNE. The data observations are divided into ten classes, each representing one digit. We sample the full MNIST dataset by selecting only 1500 random samples from each class. The projections are visualized in Fig. 8b. The PCA projection introduces significant amount of overlap between classes and the cost associated with this projection is L = 4.4 . The cost for a much more unique t-SNE map, with clearly separated classes, is lower, L = 1.7 . Figure 8c additionally shows a comparison of the D (σ ) curves corresponding to the 2D PCA projection and the 2D t-SNE map, showing that the D (σ ) metric detected the significant class overlap in PCA. With some amount of scattered www.nature.com/scientificreports/ observations still present in the t-SNE map, the D (σ ) is not composed of a single rise as was seen in the fully non-overlapping classes example in Fig. 8a. This behavior in D (σ ) can be exploited to inform the appropriate hyper-parameter tuning to improve t-SNE maps. In the Supplementary material, we further demonstrate the potential of the cost function to guide the selection of an important hyper-parameter of t-SNE called perplexity.

Discussion
Many factors can affect the quality of low-dimensional data parameterizations and there is a need to quantitatively assess those factors from the reduced-order modeling perspective. We propose a cost function that reduces the low-dimensional manifold topology to a single number. The two topological properties that the cost function pays attention to are uniqueness and feature sizes in relevant dependent variable(s). There are two main strengths of the proposed cost function. First, the manifold topology can be optimized for any target dimensionality. Second, the manifold topology can be optimized with respect to any user-selected dependent variables. Only the most important dependent variables need to be included in the manifold topology optimization. This can become particularly helpful in approaches where large state-spaces are compressed to a smaller number of parameters and regression is used to predict a physical quantity based on the compressed representation. Optimal manifolds can be found specifically from the perspective of the physical quantities of interest. We demonstrate applications on numerically generated datasets and on experimental data containing noise. Our cost function can have useful applications in searching for the best data preprocessing strategies. This can include benchmarking existing strategies with newly invented ones. Often, those settings need to be tailored to a specific dataset and even to the target manifold dimensionality. In addition, various dimensionality reduction and manifold learning techniques can be assessed from the point of view of generating quality manifolds. A possible application might be to optimize the hyper-parameters to result in improved manifold topologies, especially for techniques such as t-SNE or UMAP which are known to strongly depend on hyper-parameter settings. While visual inspection of low-dimensional manifolds can be helpful, our cost function is a quantitative metric that can aid in qualitative assessments, especially when the manifold dimensionality exceeds 3D. Good quality of manifold topology is crucial when nonlinear regression is employed on manifolds. We show that improved projections can bring modeling benefits in regression tasks. We also briefly delineate possible applications of the cost function in dealing with categorical data where the cost function is computed from a dependent variable containing discrete values (class labels). For categorical data, the numerical values for the class labels and the distances between these values will affect the cost. Future work can investigate further the application to categorical data, especially for datasets where the number of classes becomes large or where classes exhibit a meaningful hierarchical structure, such as genomics or transcriptomics data. Future work can also include incorporating the cost function directly as an objective function in dimensionality reduction.
There remain some limitations of our proposed cost function. First, the computational cost for the normalized variance derivative, D (σ ) , becomes high for large datasets. A potential solution to this restriction is to sample the dataset prior to computing the cost function. The Supplementary material sheds some light on how the cost function might respond to data sampling. Second, the current definition of the cost function cannot distinguish between the multiple scales of variation on a unique manifold and the non-uniqueness alone. This is due to the fact that the additional peaks in D (σ ) will show up in both scenarios, and both will increase the area under the D (σ ) curve. In our experience, however, the raise in D (σ ) due to non-uniqueness usually happens for σ ≪ σ peak , while the raise in D (σ ) due to varying feature sizes on an otherwise unique manifold happens at length scales closer σ peak (cf. Fig. 3d,f). The work in Ref. 30 demonstrates how peak locations in D (σ ) due to non-uniqueness show a sensitivity to data sampling that peak locations due to unique features do not. Future work can consider incorporating information such as this into the cost function to better distinguish manifolds with non-uniqueness from those with small features. Finally, the resulting cost for a given manifold can only be interpreted in relation to cost(s) for other manifold(s). Thus, the cost function can help identify the best manifold among a set of manifolds, but, as of yet, no objective judgment can be made for a single cost value obtained for www.nature.com/scientificreports/ any one manifold. The same applies to cost associated with any one dependent variable on a single manifold-it should be interpreted in relation to costs for other dependent variables. Nevertheless, we believe that this is a timely contribution that can help researchers across various disciplines and across numerous applications where low-dimensional data projections play an important role. Although we focused the demonstrations in this paper to four main datasets (argon plasma, reacting flows, atmospheric pollutant dispersion and categorical data), the cost function proposed can have broad applications in virtually any domain of science. We particularly look forward to exploring other applications, for instance on biological or medical data. We argue that further improvements in parameterization quality can be achieved in many areas of research if the low-dimensional parameter space is thoroughly explored and then assessed using the proposed cost function.

Methods
Data normalization. Prior to applying dimensionality reduction, datasets are often normalized. This can be especially beneficial if the dataset is composed of variables that have very different numerical ranges. Given a dataset composed of Q variables, X = [X 1 , X 2 , . . . , X Q ] , we normalize each variable X j by subtracting its center, c j , and dividing it by the scaling factor, d j . In a matrix form, we can write the normalized dataset, X , as: where C ∈ R N×Q is a matrix of centers with the jth column of that matrix populated with a value c j and D ∈ R Q×Q is a diagonal matrix of scales with the jth element on the diagonal equal to d j . In this work, we adopt several data scaling criteria collected in Table 1, where s j is the standard deviation and k j is the kurtosis of X j . Techniques S1-S3 are newly introduced scaling techniques that we explore in this work. Throughout the main text, we refer to a particular scaling by its name as per Table 1. The name "None" is equivalent to no scaling applied to the dataset.
Outlier removal. For the atmospheric pollutant dispersion dataset, we perform outlier detection and removal using the principal component classifier method 53,75 . Outliers are detected based on major and minor principal components (PCs). The observation i is classified as an outlier if the first PC classifier based on the q first (major) PCs: or if the second PC classifier based on the (Q − k + 1) last (minor) PCs: where z ij is the (i, j)th element from the PCA-transformed data matrix and L j is the jth eigenvalue from PCA. Major PCs are selected such that the total variance explained is 50% (this determines the number q). Minor PCs are selected such that the remaining variance they explain is 20% (this determines the number k). PCA is performed with Auto scaling. Coefficients c 1 and c 2 are found such that they represent the quantile of the empirical distributions of the first and second PC classifier respectively. Here, we set the quantile to 98%, which allowed to find 509 outlying observations out of 18,540 total observations. Dimensionality reduction. A linear projection of a dataset onto a basis defined by A ∈ R Q×q can be performed as: where η η η = [η 1 , η 2 , . . . , η q ] defines the q-dimensional manifold parameters. In PCA, which is frequently used in this work, the basis matrix A is computed as the eigenvectors of a data covariance matrix. In the equation above, we assume that the dataset, X , has already been appropriately preprocessed. The summary of linear and nonlinear dimensionality reduction techniques explored in this work is given in Table 2. For results reproducibility, random seed of 100 is used in all techniques that rely on randomness.
PCA. We use the PCA implementation from the PCAfold Python package 77 developed by the authors.

LDA.
We use the LDA implementation from the scikit-learn library. For the swiss roll dataset example, all parameters are set to default.

DML.
We use the the DML implementation from the pyDML Python package 76 . For the swiss roll dataset example, all parameters are set to default.

MDS.
We use the MDS implementation from the scikit-learn library. For the swiss roll dataset example, all parameters are set to default. Autoencoder. We use the Keras Python library to set up the autoencoder for the swiss roll dataset example.
For results reproducibility, random seed of 100 is used in TensorFlow (tf.random.set_seed(100)). The model architecture is 3-2-3 with hyperbolic tangent activations in all layers. The network weights are initialized from the Glorot uniform distribution and the biases are initially set to zeros. We use batch size of 50 and 200 epochs. The Adam optimizer is used with the learning rate 0.001 and the loss function is the mean squared error.
We train the autoencoder model on 80% of the data.
The proposed cost function. The starting point for formulating our cost function is computing the normalized variance proposed by Armstrong and Sutherland 30 . The goal is to assess the low-dimensional parameterization quality defined by q manifold parameters, η η η ∈ R N×q , obtained using any dimensionality reduction or manifold learning technique. For length scales on a manifold given by the parameter σ ∈ �σ min , σ max � , the normalized variance, N (σ ) , is computed for the ith dependent variable, φ i , as: where N is the number of observations in a dataset, φ i is an arithmetic average of φ i , and K i is computed as a weighted average of observations of φ i : where the weights, w j , are determined using a Gaussian kernel: Table 2. Dimensionality reduction and manifold learning techniques used in this work. In the last column, we list the most important hyper-parameters of each method. The cost function can help fine-tune the hyperparameters to achieve quality manifold topologies. www.nature.com/scientificreports/ In Eq. (7), the quantity ||η j − η|| 2 2 is the squared Euclidean distance between the current location on a manifold, η , and any jth point on a manifold, η j . We then construct the normalized variance derivative function as per 30 : and normalize it by its maximum value: The basis for the cost function proposed in this work is integration of the penalized normalized variance derivative function, D (σ ) , over length scales given by σ . For the ith dependent variable, φ i , the area under the D i (σ ) curve can be computed as: where the tilde denotes a log 10 -transformed quantity (e.g. σ = log 10 σ ). We compute the area in the log 10 -space of the length scales σ so that all scales of variation with different orders of magnitude are treated equally. We further introduce two penalties when computing the area: 1. Penalty for peak locations relative to the rightmost peak, σ peak . This favors large feature sizes on a manifold over small ones. 2. Penalty for the area under D (σ ) happening at σ < σ peak , and especially for σ ≪ σ peak . This penalizes multiple scales of variation that might show up as additional peaks in D (σ ) and becomes particularly useful when these additional peaks can be linked to non-uniqueness.
The cost function proposed herein takes these two penalties into account and is defined for φ i as: where P i (σ , σ peak,i ) is the penalty function defined as: The first term in P i penalizes non-zero values in D i (σ ) at length scales σ < σ peak . It will especially amplify any area under the D i (σ ) curve happening for σ ≪ σ peak . By increasing the power, the user can increase the amount of penalty for the variance occurring at σ ≪ σ peak . The second term in P i introduces a gentle penalty for the location of the rightmost peak and essentially increases or decreases the entire penalty function by a constant value. This second term acts to reward larger feature sizes as those should be easier to model. The parameter b is another hyper-parameter which controls the amount of the constant vertical shift of the entire penalty function. By increasing b, the user can increase the amount of penalty for the rightmost peak location, σ peak . Throughout this work, we use r = 1 and b = 1 . For completeness, in the next section, we illustrate the effect of setting r and b to values other than unity. Python implementation of N (σ ) , D (σ ) and the proposed cost function, L , has been developed by the authors and is available in the PCAfold software package 77 . In practice, numeric integration of D (σ ) is performed using a composite trapezoid rule. The peak values in D (σ ) are computed using the scipy.find_peaks function, which finds all local maxima by comparing the neighborhood of all discrete values of D (σ ) . The value for σ peak is taken as the rightmost peak found by scipy.find_peaks, but the user can also select a percentage, p, of the value σ peak (such that the true rightmost peak becomes p(σ max − σ peak ) ). Our recommended range for the parameter σ is �10 −7 ; 10 3 � , with logarithmically-spaced in-between values. For instance, throughout this work we typically set sigma = numpy.logspace (-7, 3, 200). Throughout this work, the manifold parameters η η η are scaled to a unit box (each η i is scaled to a [0, 1] range) before computing the cost. This is done so that the length scale, σ , has the same meaning in each manifold dimension.
Effect of hyper-parameters r and b on the cost function. The hyper-parameters r and b from Eq. (12) may be changed to emphasize penalties applied to non-uniqueness and/or small feature sizes when computing overall costs. Figure 9 illustrates the effect of setting r and b to non-unity values. In Fig. 9a-c, we utilize the toy functions from Fig. 3 and show how L behaves for various r and b. For clearer analysis, whenever r is varied, b is set to unity, and vice versa. For comparison, with the red dashed lines, we mark costs corresponding to r = 1 and b = 1 . In principle, increasing b increases L more severely for manifolds where a dependent variable exhibits small feature sizes (see Fig. 9a). Conversely, increasing r for unique manifolds does not affect the L values significantly. As can be understood from Fig. 9c, increasing r increases L more severely for manifolds with ) . www.nature.com/scientificreports/ non-uniqueness. Even if non-uniqueness is present on a manifold with a fixed largest feature size, changing b does not affect the L values significantly. For a unique manifold with multiple feature sizes, increasing r or b can increase L by a similar amount (Fig. 9b). Due to the blending of multiple scales of variation in the D (σ ) profile seen in Fig. 3d, we have also applied a p = 70% shift in σ peak in Fig. 9b, as per discussion in the previous section. In Fig. 9a-c, the overall trend of the cost function is preserved when changing r and b. We further test the impact of r and b on the verdict given by the cost function across various dimensionality reduction and manifold learning techniques. We use the classic 3D swiss roll dataset and generate various 2D projections. The original 3D topology of this dataset can be seen in Fig. 9d, colored by a dependent variable, φ . The cost corresponding to the parameterization in the original 3D space is L = 0.98 with r = 1 and b = 1 . Figure 9e shows twelve different 2D projections of that dataset generated with various linear and nonlinear techniques. The cost value, L , is reported for each projection taking r = 1 and b = 1 . Among all the reduction techniques selected for this demonstration, PCA and a linear autoencoder (AE) introduce the most significant non-uniqueness. The costs for AE and PCA projections are the highest, while SE, LDA and UMAP show the lowest costs. The small cost for SE, LDA or UMAP projections (comparable to the cost in the original 3D space) is likely due to a combination of two factors: projection uniqueness and large feature sizes created through sufficient separation of distinct φ values in space. Figure 9f shows the D (σ ) curves for the original 3D data (red dashed lines) and for the 2D projections, in groups of four. The location of σ peak is shifted to the right for the SE and UMAP projections with respect to the original 3D data parameters. A similar observation holds for LLE, H-LLE, LTSA and Isomap projections. This can be due to a relatively large distance between the smallest and the largest φ values on these projections. In the original 3D space, the high and low values of φ are relatively close to each other due to manifold curvature. In Fig. 9g,h, we show the effect on the cost function evaluations of the various swiss roll data projections when changing the hyper-parameters r and b. In Fig. 9g, costs increase with increasing r only for the PCA and AE projections. This is understandable, since those are the only projections with significant levels of non-uniqueness. Increasing the hyper-parameter r thus increasingly amplifies the "secondary" rise seen in the D (σ ) profile at σ < σ peak for PCA and AE. Interestingly, this behavior with r can potentially be used to detect manifolds with non-uniqueness among a set of parameterizations. In Fig. 9h,  www.nature.com/scientificreports/ increasing b increases costs for all parameterizations (3D and 2D). Costs increase more rapidly for PCA and AE projections than for any other explored projection (note the logarithmic scale of the vertical axis). This is also due to emphasizing the "secondary" rise seen in the D (σ ) profile when computing the penalized area, but this time through increasing the vertical shift of the entire penalty function. With the circled outlines in Fig. 9g,h, we mark the lowest cost that happened for any of the 2D projections across the explored values of r and b. Among the r and b values explored, the smallest cost consistently happens for the SE projection. We also note that the ratio between L for the worst projection (AE) and L for the best projection (SE) is generally higher when r is set above unity than when b is set above unity. Thus, increasing r above unity can help create clearer separations in L values when ranking manifolds with varying levels of non-uniqueness.

(10)
On the computational impact of evaluating the cost function. We note that the computational time required to compute L for a dataset with N observations and for m dependent variables scales as O(mN 2 ) , assuming the same number of discrete values of σ . The code for computing the cost function, that we provide in the PCAfold library, is parallelized with respect to σ , since the normalized variance computations are entirely independent of one another for different values of σ . Thus, our code can be readily ran on multiple CPUs, which we generally recommend for N > 10 5 . In the Supplementary material, we show the effect of data sampling on the cost function. We note, that while subsampling the data can ease the computational cost of L , it has an enhanced effect on dependent variables which exhibit variation at multiple scales (e.g. from non-uniqueness).
A possible future implementation that could reduce the computational time of evaluating Eq. (6) is to approximate the weighted average with a sum over points within/near the currently considered length scale, σ . This implementation can make the computational time scale approximately as O(mN) for small σ and only scale as O(mN 2 ) as σ encompasses all points on a manifold. There is a second aspect to the question of the computational cost, that is more elusive to quantify, and that is the time required for the researcher to find a good manifold topology through trial and error exploration, if no quantitative tools were applied to guide the choice of a manifold. We argue that this second aspect can become a bottleneck in effective application of reduced-order models. In this regard, one may find that our cost function allows for "quick" assessments of projections resulting from a range of dimensionality reduction and manifold learning techniques, as well as a variety of data preprocessing strategies applied to the training data. Moreover, since L can be implemented in optimization tasks, it is likely that an optimum in manifold topologies is found in an automated way, compared to a trial and error approach.
L-informed feature selection. We developed a feature selection algorithm that iteratively eliminates state variables from the dataset based on minimizing the cost, L , of PCA projections 65 . There are three inputs to the algorithm: the original dataset, X ∈ R N×Q , the target dependent variables, φ φ φ ∈ R N×m , and the target manifold dimensionality, q. The pseudocode below shows in closer detail how the algorithm computes the optimized subset of the original state variables. At each iteration i, the algorithm computes PCA projections resulting from removing each variable, one at a time. At the end of the iteration, the variable whose removal decreased the cost value the most is discarded from the dataset. At the next iteration, the process repeats, but now on a dataset with one less variable. We only allow Q − q iterations so that we never reduce original data dimensionality below q requested by the user. Once all iterations have finished, the algorithm looks back at final costs from all iterations and returns the optimized subset, X S , corresponding to the iteration that showed the minimum cost value. Such subset should then generate an optimized manifold topology. The algorithm is available in the PCAfold software package 77 .
Nonlinear regression using artificial neural networks (ANNs). In the atmospheric pollutant dispersion example, we use nonlinear regression using ANNs. We use the Keras Python library to set up the ANN regression. For results reproducibility, random seed of 100 is used in TensorFlow (tf.random.set_ seed(100)). The inputs of the network are the three manifold parameters, η η η = [η 1 , η 2 , η 3 ] , either from PCA, UMAP, Isomap or t-SNE. The output is the turbulent Schmidt number, φ = Sc t . The model architecture is 2-5-5-1 for 2D projections and 3-5-5-1 for 3D projections with sigmoid activations in all layers except the output layer, where we use linear activation. The network weights are initialized from the Glorot uniform distribution www.nature.com/scientificreports/ and the biases are initially set to zeros. We use batch size of 100, validation split of 0.2 and 500 epochs. The Adam optimizer is used with the learning rate 0.001 and the loss function is the mean squared error. We train the regression model on 80% of the data and measure the MAE on the remaining 20% test data not seen by the ANN model. The results reported in Fig. 7 are for the 20% test data.
Nonlinear regression using kernel regression. In the atmospheric pollutant dispersion example, we use kernel regression from the PCAfold software package 77 . We use the Nadaraya-Watson estimator with a Gaussian kernel with a varying bandwidth. The bandwidth is determined locally based on 50 nearest neighbors of the query point. We train the regression model on 80% of the data and measure the MAE on the remaining 20% test data not seen by the kernel regression model. The results reported in Fig. 7 are for the 20% test data.
Reacting flow data generation. The reacting flow datasets for combustion of syngas in air and hydrogen in air were generated using Spitfire Python package 78 available at: github.com/sandialabs/Spitfire. The datasets were generated using a steady laminar flamelet model for a range of dissipation rates from chemical equilibrium to extinction and a range of mixture fractions between 0 and 1. For the syngas/air dataset, the fuel stream is composed of a mixture of carbon monoxide and hydrogen in 10:1 molar proportion 79 . The oxidizer stream is air. Both streams have initial temperature 300 K and pressure 101,325 Pa. For the hydrogen/air dataset, the fuel stream is composed of hydrogen and the oxidizer stream is air 80 . Both streams have initial temperature 300 K and pressure 101,325 Pa. For both fuels, the mixture fraction grid is denser closer to the stoichiometric conditions.